Integration of visual information, anatomic constraints and prior shape knowledge for medical segmentations

ABSTRACT

This invention relates to the integration of visual information, anatomic constraints, prior shape knowledge, and level set representations for the segmentation of medical images. An embodiment according to the present invention comprises a level set variational framework that uses a bi-directional boundary flow, an intensity-based regional component that maximizes the a posteriori segmentation probability, a physiology-based module that constrains the solution space and a term that accounts for shape-driven consistency. All modules are expressed in an energetic form and the resulting objective function is optimized using a gradient descent method.

CROSS REFERENCE TO RELATED APPLICATION

This application claims the benefit of U.S. Provisional Application Ser.No. 60/354,005 filed on Feb. 1, 2002, and U.S. Provisional ApplicationSer. No. 60/354,004 filed on Feb. 1, 2002, which are incorporated byreference herein in their entirety.

FIELD OF THE INVENTION

This invention relates to the integration of visual information,anatomic constraints, prior shape knowledge, within level setrepresentations for medical segmentation and grouping.

BACKGROUND OF THE INVENTION

Cardiovascular diseases are one of the leading causes of death in UnitedStates. The development of new technologies to diagnose cardiovasculardiseases has led to a decline of the mortality rate. Magnetic resonanceimaging (“MRI”) provides time varying three dimensional imagery of theheart that can be processed using Computer Vision techniques. The threedimensional imagery can then be used for diagnostic purposes. Therefore,physicians are interested in identifying heart chambers using threedimensional imagery so that this information can be used for diagnosticpurposes.

In medical applications the areas of interest in the image domain areactual anatomic structures. The segmentation result is constrainedaccording to some a priori high level knowledge of an anatomicstructure, such as, forms, relative positions, and motion over time.

The extraction of accurate results is a requirement in medicalapplications. Therefore, predominately boundary-based methods, enforcedby some region-based segmentation modules, have been employed to copewith the segmentation task. These methods can be classified into twocategories. The first category is parametric methods that determine thesegmentation map by fitting a boundary template to the image. The secondcategory is non-parametric methods that are based on the propagation ofregular curves under the influence of local image characteristics.

Parametric methods have real-time performance and can deal withincomplete data and slight deformations due to shape-driven constraints.However, such methods primarily refer to boundary information andrequire complicated models to deal with important shape deformationswhile topological changes cannot be handled.

Non-Parametric methods can deal with important shape deformations andtopological changes. At the same time these methods do not require apriori knowledge. However, non-parametric methods do not have a robustbehavior due to the presence of noise and are time consuming.

Identifying heart chambers, in particular, the endocardium and theepicardium, is a challenging problem in medical image analysis.Furthermore, measuring ventricular blood volume, ventricular wall mass,ventricular wall motion and wall thickening properties over variousstages of a cardiac cycle are powerful diagnostic tools that are studiedin medical image analysis. The left ventricle is of particular interestbecause it pumps oxygenated blood out of the heart to distant tissue inthe entire body. The information space refers to physically corrupteddata. Small parts, for example, papillary muscles, of an endocardium donot have the same intensity properties as its dominant parts. At thesame time, the separation, that is, boundaries, of the epicardium andthe rest of the heart is practically impossible to detect in some areas,due to lack of information.

A need exists for a unified mathematical framework that makes optimaluse of visual information, translates high level application constraintsinto low level segmentation modules, is able to deal with lack of visualinformation, and is able to deal with physically corrupted data.

SUMMARY OF THE INVENTION

An embodiment according to the present invention allows for identifyingheart chambers, in particular, the endocardium and the epicardium.Another embodiment according to the present invention measuresventricular blood volume, ventricular wall mass, ventricular wall motionand wall thickening properties over various stages of a cardiac cycle.

An embodiment according to the present invention further comprises twoprinciples. The first principle is that a segmentation map is inaccordance with visual information, such as, boundary and regionalinformation. The second principle is that the segmentation map has torespect the physiology of anatomical structures, that is, relativepositions and shape forms, of the heart. Therefore, an embodimentaccording to the present invention comprises a level set variationalframework that uses a bi-directional boundary flow, an intensity-basedregional component that maximizes the a posteriori segmentationprobability, a physiology-based module that constrains the solutionspace and a term that accounts for shape-driven consistency. All modulesare expressed in an energetic form and the resulting objective functionis optimized using a gradient descent method.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow diagram illustrating a method using anatomy shapepriors for level set representation according to an embodiment of thepresent invention.

FIG. 2 illustrates is a functional block diagram according to anembodiment of the present invention.

FIG. 3 is an image illustrating an epicardium interface and anendocardium interface according to an embodiment of the presentinvention.

FIG. 4 is an image illustrating MR sequences of a left ventricleaccording to an embodiment of the present invention.

FIG. 5 illustrates a coupling function according to an embodiment of thepresent invention.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

Segmenting medical structures is a challenging application in computervision. Referring to FIG. 1, a method according to an embodiment of thepresent invention comprises an energetic form,that introduces shapeconstraints to level set representations. Step 102 includes developing amedical shape model on level set space for tracking moving interfaces.The development of the medical shape model is described in detail belowin Medical Shape Prior Model Construction. In step 104, visualinformation is analyzed and boundary information is integrated with themedical shape model. The integration of boundary information isdescribed in detail below in Boundary Module. In step 106, visualinformation is analyzed and regional information is integrated with themedical shape model. The integration of regional information isdescribed in detail below in Regional Intensity-Based Module. In step108, a further integration of anatomical constraints occurs with themedical shape model. The integration of anatomical constraints isdescribed in detail below in Anatomical Constraints. In step 110, themedical shape model is integrated with physiological shape priors torecover an object of interest. The integration of the medical shapemodel with physiological shape priors is described below in Level SetPhysiological Shape Priors. This formulation of steps 102 through 110exploits all advantages of level set representations resulting in amethod that can deal with a large number of parametric, as well as,continuous transformations. Furthermore, the integration of existinglevel set based segmentation methods leads to paradigms that can dealwith noisy, occluded, missing and physically corrupted data.

Referring to FIG. 2 a process is illustrated according to an embodimentof the present invention. A computer system 202 having a hardwarecomponent 204, for example, a server and storage device, and softwarecomponent 206, for example, an operating system and computer programs,according to an embodiment of the present invention receives input ofmathematical functionals 208, physiological shape priors 210, and amedical image 212 that is to be recovered. Shape prior propagation andminimization of non-stationary distance are used to produce recoveredmedical image 214. The process is now described in detail.

An embodiment according to the present invention comprises a variationalformulation based on level set representations to deal with thesegmentation of the left ventricle in Magnetic Resonance Images (“MRI”).The variational formulation integrates visual information, anatomicalconstraints and prior knowledge of cardiac shapes. The visualinformation is expressed through a boundary attraction component and anintensity-based grouping module. The anatomical constraints are derivedfrom physiology and refer to relative positions of cardiac structures ofinterest. Finally, shape knowledge is obtained through a construction ofa cardiac prior model that constrains the solution space within a familyof shapes. The resulting paradigm can accommodate noisy, incomplete,physically corrupted data, and cardiac deformations, such as systole anddiastole.

In an embodiment according to the present invention, extraction of themyocardium, can be achieved using three modal image segmentation.Segmentation refers to the endocardium, the epicardium and background.In particular, the area between the endocardium and the epicardium isevaluated. Using a common technique in Computer Vision, such as,evolving interfaces according to some data driven information area canbe recovered. An embodiment of the present invention can best beintroduced using the following definitions, referring to FIG. 3, ∂R_(O)represents the epicardium interface (light gray contour) 302 and ∂R₁represents the endocardium interface (dark gray contour) 304. Thesecontours 302, 304 can be used to define a three modal image partitiongiven by (i) R_(I) endocardium, (ii) R_(O)−R_(I) myocardium, and (iii)Ω−R_(O) background. The region of interest is the myocardium.

Level Set Representations

Level Set Representations are a powerful mathematical tool able to dealwith applications that share a common concern, such as, evolvinginterfaces. These representations are best understood by using a twodimensional case and an evolving interface [∂R] in an Euclidean plane,such that [∂R: [0,1]→R², p→∂R(p)]. Presuming the propagation of ∂R(p,t)takes place in the direction of its inward normal N according to ascalar function of curvature K, then

$\frac{\mathbb{d}\;}{\mathbb{d}t}\partial\left( {{p\;\_} = {{F\left( {K(p)} \right)}{N(p)}}} \right.$Then, Φ: Ω→R⁺ is a Lipchitz function that refers to a level setrepresentation:

${\Phi\left( {x,{y;t}} \right)} = \left\{ \begin{matrix}{0,{\left( {x,y} \right) \in {\partial{R(t)}}}} \\{{{+ {D\left( {\left( {x,y} \right),{\partial(t)}} \right)}} > 0},{\left( {x,y} \right) \in \;{R(t)}}} \\{{{- {D\left( {\left( {x,y} \right),{\partial(t)}} \right)}} < 0},{\left( {x,y} \right) \in \;\left\lbrack {\Omega - {R(t)}} \right\rbrack}}\end{matrix} \right.$where an evolving interface [∂R] is represented as the zero iso-surfaceof Φ and D((x,y), ∂R(t)) is the minimum Euclidean distance between thepixel (x,y) and the interface R(t) at time t.

This representation evolves in the following manner:

${\frac{\mathbb{d}\;}{\mathbb{d}t}{\Phi(p)}} = {{- {F\left( {K(p)} \right)}}{{\nabla{\Phi(p)}}}}$

These representations are implicit, parameter free, topology free andprovide a natural way to estimate geometrical properties of the evolvinginterface.

Two level set representations have to be evaluated for the segmentationof the Left Ventricle, one for the endocardium [Φ_(I)] and one for theepicardium [Φ_(O)]. The notation can be simplified as:

$\Phi_{I} = \left\{ {{\begin{matrix}{0,{\partial R_{I}}} \\{{{+ {D\left( {\partial R_{I}} \right)}} > 0},R_{I},} \\{{{- {D\left( {\partial R_{I}} \right)}} < 0},\left\lbrack {\Omega - R_{I}} \right\rbrack}\end{matrix}\mspace{20mu}{\Phi_{O}(p)}} = \left\{ \begin{matrix}{0,{\partial R_{I}}} \\{{{+ {D\left( {\partial R_{o}} \right)}} > 0},R_{o},} \\{{{- {D\left( {\partial R_{o}} \right)}} < 0},\left\lbrack {\Omega - R_{o}} \right\rbrack}\end{matrix} \right.} \right.$

To connect the use of these representations with energy minimizationtechniques and the propagation of curves, define the approximation ofDirac and Heaviside distributions as:

${\delta_{\alpha}(\phi)} = \left\{ {{\begin{matrix}{0,{{\phi } > \alpha}} \\{{\frac{1}{2\alpha}\left( {1 + {c\; p\;{s\left( \frac{\pi\;\phi}{a} \right)}}} \right)},{{\phi } < \alpha}}\end{matrix}{H_{\alpha}(\phi)}} = \left\{ \begin{matrix}{1,{\phi > \alpha}} \\{0,{\phi < {- \alpha}}} \\{{\frac{1}{2}\left( {1 + \frac{\phi}{\alpha} + {\frac{1}{\pi}{\sin\left( \frac{\pi\;\phi}{\alpha} \right)}}} \right)},{{\phi } < \alpha}}\end{matrix} \right.} \right.$

Then it can be demonstrated that{(x,y)εΩ: lim_(a→o) ₊ [H _(a)(φ((x,y);t))]=1}=R{(x,y)εΩ: lim_(a→o) ₊ [δ_(a)(φ((x,y);t))]=1}=∂R.

Using these components, define boundary-based modules and region-basedmodules for the evolving interface as follows:

$\begin{matrix}{\underset{\underset{{region}\mspace{14mu}{module}}{︸}}{\int{\int_{\Omega}^{\;}{{H_{\alpha}\left( {\Phi\left( {x,y} \right)} \right)}{r\left( {I\left( {x,y} \right)} \right)}\ {\mathbb{d}x}{\mathbb{d}y}}}}\mspace{95mu}} & (i) \\{\underset{\underset{{boundary}\mspace{14mu}{module}}{︸}}{\left. {\int{\int_{\Omega}^{\;}{{\delta_{\alpha}\left( {\Phi\left( {x,y} \right)} \right)}{b\left( {I\left( {x,y} \right)} \right)}}}} \middle| {{\nabla{\Phi\left( {x,y} \right)}}\ {\mathbb{d}x}{\mathbb{d}y}} \right.}\;} & ({ii})\end{matrix}$where r and g are region and boundary attraction functions. Thesemodules have the following interpretation: the first term [i] refers toa grouping component since it accounts for some regional properties,that is modulo the definition of r, of the area defined by the evolvinginterface. The second term [ii] refers to the interface points and theterm can be a combination of a boundary attraction term, that is, modulothe definition of b, and a smoothness component, such as, a geodesicactive contour.

These general purpose segmentation modules can now be used to provide adata-driven solution for the extraction of the myocardium.

Visual Information Modules

Image segmentation approaches can be classified as eitherboundary-based, that rely on the generation of a strength image and theextraction of prominent edges, or region-based, that rely on thehomogeneity of spatially localized features and properties.

Boundary Module

A common way to determine boundary information is through a gradientimage and a monotonically decreasing function g( ), for example,

${{g(I)} = \frac{1}{1 + {{\nabla I_{s}}}}},{{g(I)} = {\frac{1}{\sqrt{2\;\pi}\sigma}{{\mathbb{e}}^{- \frac{{\nabla I_{s}}}{2\sigma^{2}}}.}}}$Then, within the level set representation, a geodesic active contourmodel can be considered to provide a boundary-based segmentationsolution. The geodesic active contour model is directed to finding aminimal length geodesic curve taking into account desired imagecharacteristics. Given the application framework, two energycomponents/geodesic curves, such as endocardium and epicardium, areused:

${E\left( {\Phi_{I},\Phi_{O}} \right)} = {\underset{\underset{endocardiumboundaryattraction}{︸}}{\int{\int_{\Omega}^{\;}{{\delta_{\alpha}\left( \Phi_{I} \right)}{g(I)}{{\nabla\Phi_{I}}}}}} + \underset{\underset{epicardiumboundaryattraction}{︸}}{\int{\int_{\Omega}^{\;}{{\delta_{\alpha}\left( \Phi_{O} \right)}{g(I)}{{\nabla\Phi_{O}}}}}}}$

This energy is minimized using a gradient descent method. The calculusof variations provides flows that shrink initial interfaces towards thecardiac boundaries while respecting some internal regularityconstraints. Although these flows have been proven to be efficient innumerous applications, they suffer from being myopic since the datadriven term tends to shrink the interface, such as in single directionalflows. Although this component may not be appropriate for boundaryextraction, it can be used to introduce constraints that account forsome desired internal properties, for example, smoothness, of theevolving interfaces:

${E\left( {\Phi_{I},\Phi_{O}} \right)} = {\underset{\underset{endocardiumsmoothness}{︸}}{\int{\int_{\Omega}^{\;}{{\delta_{\alpha}\left( \Phi_{I} \right)}{{\nabla\Phi_{I}}}}}} + {\underset{\underset{epicardiumsmoothness}{︸}}{\int{\int_{\Omega}^{\;}{{\delta_{\alpha}\left( \Phi_{O} \right)}{{\nabla\Phi_{O}}}}}}.}}$

The limitation of the Geodesic Active Contour model, that is, a singledirectional flow, can be dealt with using gradient vector flow (“GVF”).This field refers to the definition of a bi-directional external forcethat captures object boundaries from either side and can accommodateconcave regions. The GVF comprises a two dimensional vector field[v(p)=(V^(x)(p), V^(y)(p)), p=(x,y)] that minimizes the followingenergy:

E(v) = ∫∫μ((v_(x)^(x))² + (v_(y)^(x))² + (v_(x)^(y))² + (v_(x)^(y))²) + ∇f²v − ∇f²𝕕x𝕕y,where [ƒ(p)=1−g(p)] and μ is a blending parameter. According to thisfunction, areas where the information is constant [|∇ƒ|→0], aredominated by partial derivatives of the field, such as, smooth flow map.Discontinuities of the boundary space, that is |∇ƒ| is large, activatethe second term energy, leading to v=∇ƒ.

The objective function does not make direct use of the boundaryinformation, that is, only the gradient of the boundary affects theflow. This characteristic can be considered a limitation since strongedges, as well as, weak edges may create a similar flow due to thediffusion process. To overcome this problem, the objective function canbe modified by introducing boundary information in a direct form, suchas:E(v)=∫∫μ(v _(x) ^(x))²+(v _(y) ^(x))²+(v _(x) ^(y))²+(v _(x)^(y))²)+ƒ|∇ƒ|² |v−∇ƒ| ² dxdyThis modification induces the ability of overcoming weak edges due tonoise presence. Also, it leads to an appropriate diffusion of theboundary information. Strong edges overcome/compensate flows that areproduced by weak edges. The minimization of the objective function willlead to a diffusion flow that propagates the information from theobjects boundaries to the background.

This field was originally used to propose a bi-directional flow that canreach an object's boundaries from either side. However, changingtopology was an important limitation of the original flow. Changingtopology was dealt with in a non-energetic form by integrating the GVFto the geodesic active contour flow.

An embodiment according to the present invention integrates the GVF flowto level set representations in an energetic form. The increase of fluxof an auxiliary vector field can be optimized using images of lowcontrast, for example, elongated structures, such as blood vessels.Thus, the following minimization criterion can be achieved:

$E_{B}\left( {\Phi_{I},{\Phi_{O} = {\underset{\underset{endocardiumboundaryattraction}{︸}}{\int{\int_{\Omega}{{\delta_{\alpha}\left( \Phi_{I} \right)}\left\lbrack {{\nabla\Phi_{I}} \cdot \left( {v^{x},y^{y}} \right)} \right\rbrack}}} + \underset{\underset{epicardiumboundaryattraction}{︸}}{\int{\int_{\Omega}{{\delta_{\alpha}\left( \Phi_{O} \right)}\left\lbrack {{\nabla\Phi_{O}} \cdot \left( {v^{x},y^{y}} \right)} \right\rbrack}}}}}} \right.$The minimization of this energy component can be done using a gradientdescent method and the calculus of variations. The obtained motionequations refer to a bi-directional flow that can reach the object'sboundaries from either side.Regional Intensity-Based Module

Regional/global information has been used to provide segmentationresults and improve the performance of boundary-based flows. A regionalintensity module can use an evolving interface to define an imagepartition that is optimal with respect to some grouping criterion. Formedical images, observed intensities depend on the properties of thecorresponding tissue being mapped.

Furthermore, referring to FIG. 4, the MR sequences of the left ventriclehave three populations: (i) blood (bright) 402, (ii) muscles (gray) 404,and (iii) air-filled lungs (dark gray) 406. The characteristics of thesepopulations vary spatially and temporally, but their intensityproperties can be discriminated according to their relative differences.Therefore, an observed distribution, that is, a histogram, of a wideepicardium region is a Gaussian mixture model with three components. Letp_(I) be the endocardium density function, p_(O) the myocardium densityfunction and p_(B) the density function of the rest of the cardiacorgans, that is the background.

Then,p(I)=p _(I) p _(I)(I)+p _(O) p _(O)(I)+p _(B) p _(B)(I)where p_(I), p_(O), p_(B) are the a priori probabilities for theendocardium, the myocardium, and the background components. The unknownparameters of this model can be estimated using anexpectation-maximization principle.

The image partition that accounts for the expected image characteristicsof the different components is the one that maximizes the a posteriorisegmentation probability:

${E\left( {\Phi_{I},\Phi_{O}} \right)} = {\underset{\underset{endocardium}{︸}}{\int{\int_{\Omega}{{H_{O}\left( \Phi_{I} \right)}{g\left( {p_{I}(I)} \right)}}}} + \underset{\underset{mycardium}{︸}}{\int{\int_{\Omega}{{H_{\alpha}\left( \Phi_{O} \right)}\left( {1 - {H\;{\alpha\left( \Phi_{I} \right)}}} \right){g\left( {{po}(I)} \right)}}}} + \underset{\underset{background}{︸}}{\int{\int_{\Omega}{\left( {1 - {H\;{\alpha\left( \Phi_{O} \right)}}} \right){g\left( {p_{B}(I)} \right)}}}}}$where g is a monotonically decreasing function. The myocardium is theregion between the epicardium [Φ_(O), H_(α) (Φ_(I))>0] and theendocardium [Φ₁, H_(α) (Φ_(I))<0].

Although the use of this module improves the performance of thesegmentation algorithm, there are cases where it produces sub-optimalresults. The intensity properties of the papillary muscles do not followthe endocardium distribution and most of the cardiac organs located atthe bottom-left part of the left ventricle follow the myocardiumintensity properties. Therefore, the papillary muscles cannot besegmented as part of the endocardium, while other cardiac organs areconsidered to be part of the myocardium. There is significant overlapbetween the Gaussian component of the myocardium and the background. Toadjust for these cases, anatomy-based and shape-drive constraints haveto be introduced.

Anatomy and Shape-Driven Modules

Using visual information has been proven to be very efficient in imagesegmentation/grouping applications. However, for the left ventriclesegmentation problem solving the inference based on visual informationis not always feasible. The physiology of the heart and its behaviorduring cardiac cycles, has been exhaustively studied. Therefore,knowledge from other domains have to used to derive low level modulesthat improve the segmentation performance of the algorithm.

Anatomical Constraints

According to visual information, segmentation of the left ventricle incardiac MR images has been associated with two evolving interfaces, thatis, inner cardiac contours and outer cardiac contours.

In the left ventricle case, evolving interfaces are used to describe thebehavior of some components of the heart. Their relative positions areconnected in some way and have to respect some given anatomicalconstraints.

Therefore, distance and respectively, the relative positions of theevolving interfaces, can be constrained to be within a given set ofacceptable values. Similar constraints that have been motivated by humananatomy have been used in the prior art for the segmentation of thecortex and the left ventricle. Evolving motion equations that preservethe distance constraint have been proposed and also introduced using anenergy minimization framework.

The use of level set representations and the distance transform as anembedding function results in a straight-forward way to estimate thedistance between evolving interfaces. Referring to FIG. 5, a couplingfunction h is shown as:

${h\;{\delta(x)}} = \left\{ \begin{matrix}{1,{{{if}\mspace{14mu}\left\lbrack {x \leq n} \right\rbrack}\bigcup\left\lbrack {x > M} \right\rbrack}} \\{0,{{if}\mspace{14mu}\left\lbrack {m < x < M} \right\rbrack}} \\{\left( \frac{x - m}{\delta} \right)^{2},{{if}\mspace{14mu}\left\lbrack {{m - \delta} < x \leq m} \right\rbrack}} \\{\left( \frac{x - M + \delta}{\delta} \right)^{2},{{if}\mspace{14mu}\left\lbrack {M \leq x \leq {M + \delta}} \right\rbrack}}\end{matrix} \right.$where m 502 and M 504 are the minimum and maximum allowed distancebetween the inner and the outer contour, and δ 506 is a parameter thatdetermines the extremum and the infimum of the constraint.

An anatomy-driven constraint that aims at preserving the distancebetween the endocardium and the epicardium can then be introduced. Thisconstraint also accounts for some internal properties of the evolvinginterfaces, such as, smoothness, and is given by:E_(A)(Φ_(I),Φ_(O)=∫∫_(Ω)δ_(α)(Φ_(I))h(|Φ_(O)|)|∇Φ_(I)|+∫∫_(Ω)δ_(α)(Φ_(O))h(|Φ_(I)|)|∇Φ_(O)|The interpretation of this function is clear; when the distance betweenthe evolving interfaces is within an acceptable range, this componentbecomes inactive. In contrast, when the constraint is not satisfied,then the component tends to propagate the evolving interfaces towardsthe direction that preserves the anatomical constraint. This term issimilar to the original boundary attraction component, that is, ageodesic active contour. However, it has an adaptive behavior since thedesired image properties, such as, level set representations/distancevalues, are dynamically updated.Level Set Physiological Shape Priors

Shape-driven approaches are a common selection for medical imagesegmentation. The use of B-Splines, Deformable templates, and Fourierdescriptors can provide fairly good solutions for noisy, incomplete andphysically corrupted data, such as, medical volumes. A first steptowards using shape-driven constraints refers to the definition of priorshape models, that is, parameterization of an evolving interface. Thisstep is also known as a training phase and creates a representativeshape model from a given set of examples. The present invention employsa shape-driven approach based on a shape-driven principle described in“Shape Priors for Level Set Represetantions,” by Paragyios, et al., U.S.patent application Ser. No. 60/354,005, which is incorporated byreference herein in its entirety.

Medical Shape Prior Model Construction

The next step of an embodiment according to the present invention, isconstruction of a shape model, using aligned contours. A stochasticframework with two unknown variables can be used. The variables can be ashape image, Φ_(M)(x,y) and local degrees, such as, variability, ofshape deformations σ_(M)(x,y). Each grid location can be described in ashape model using a Gaussian density function. Furthermore, a constraintthat the shape image is a level set representation can be imposed. Anembodiment according to the present invention further comprises atwo-step approach based on variational principles. The two-step approachis used to obtain the level set representation. During a first step, anoptimal solution according to data driven terms, that minimizes thedistance between the model and the input shapes, that is, maximizationof the joint density while preserving a smooth variability map, suchthat:

$\left. {{E\left( {\Phi_{M},\sigma_{M}} \right)} =} \right)\left( {1 - \alpha} \right)\begin{matrix}{\underset{\underset{smoothness}{︸}}{\int{\int_{\Omega}{\left( {\left( {\frac{\mathbb{d}\;}{\mathbb{d}x}{\sigma_{M}\left( {x,y} \right)}} \right)^{2} + \left( {\frac{\mathbb{d}\;}{\mathbb{d}y}{\sigma_{M}\left( {x,y} \right)}} \right)^{2}} \right){\mathbb{d}x}{\mathbb{d}y}}}} +} \\\underset{\underset{{data}\mspace{14mu}{attraction}}{︸}}{\frac{\alpha}{n}{\int{\int_{\Omega}{\sum\limits_{i = 1}^{n}{\left( {{\log\left\lbrack {\sigma_{M}\left( {x,y} \right)} \right\rbrack} + \frac{\left. {{\Phi_{i}\left( {x,y} \right)} - {\Phi_{M}\left( {x,y} \right)}} \right)^{2}}{2{\sigma_{M}^{2}\left( {x,y} \right)}}} \right){\mathbb{d}x}{\mathbb{d}y}}}}}}\end{matrix}$is determined. In a second step, an optimal projection of the solutionof the first step, at the manifold of acceptable solutions, such as,distance functions, is determined:

${\frac{\mathbb{d}\;}{\mathbb{d}t}\Phi_{M}} = {\left( {1 - {{sgn}\left( \Phi_{M}^{0} \right)}} \right)\left( \left. {1 -} \middle| {\nabla\Phi_{M}} \right| \right)}$The two steps alternate until the system reaches a steady-statesolution. Upon convergence of the system, a level set representation isobtained that optimally expresses the training set, using some localdegrees of variability that are spatially smooth.Level Set Shape Prior Constraint

To evolve an interface and consequently, its level set representation[Φ], while respecting some known shape properties Φ_(M)(x,y) is anobjective of a shape prior constraint. Assume that all instances of theevolving representation belong to the family shapes that is generated byapplying all possible global transformations to the prior shape model.

Thus, given Φ, there is an ideal transformation A=(A_(x),A_(y),) betweenthe shape prior and the evolving representation that maps each value tothe most probable value on the model:

$\quad\left\{ \begin{matrix}\left. \left( {x,y} \right)\rightarrow{A\left( {x,y} \right)} \right. \\{{\left. {\max_{x,y}\left\{ {p_{A}^{M}\left( {x,y} \right)} \right)} \right\}{\forall\left( {x,y} \right)}}:{{H_{\alpha}\left( {\Phi\left( {x,y} \right)} \right)} \geq 0}}\end{matrix} \right.$The most probable transformation is the one that maximizes joint densityfor all pixels. Presuming that these densities are independent acrosspixels, the minimization of the [−log] function can be considered asglobal optimization criterion:

$\left. {{E\left( {\Phi,A} \right)} = {\int{\int_{\Omega}{{H_{\alpha}(\Phi)}\left\lbrack {{\log\left( {\sigma_{M}(A)} \right)} + \frac{\left( {{s\;\Phi} - {\Phi_{M}(A)}} \right)^{2}}{2{\sigma_{M}^{2}(A)}}} \right)}}}} \right\rbrack$The interpretation of this functional is straightforward. Atransformation and a level set representation, that maximize posteriorprobability pixel-wise given the shape prior model, have to bedetermined. This model refers to a non-stationary measurement wherepixels are considered according to the confidence of their projectionsin the shape prior model.

For a left ventricle solution, presume the existence of two shape priormodels, one for the endocardium [Φ_({I,M},)σ_({I,M})] and one for theepicardium [Φ_({O,M}),σ_({O,M})]. These models can be built using manualsegmented data. Introduce a registration/shape prior component to asegmentation process that (i) registers the evolving representations tothe prior models, and (ii) deforms the evolving representations locallytowards a better matching with the prior models. Therefore,

${E_{s}\left( {\Phi_{I},\Phi_{O},A_{I},A_{O}} \right)} = {\underset{\underset{{endocardium}\mspace{14mu}{shape}\mspace{14mu}{prior}}{︸}}{\left. {\int{\int_{\Omega}{{H_{\alpha}\left( \Phi_{I} \right)}\left\lbrack {{\log\left( {\sigma_{\{{I,M}\}}\left( A_{1} \right)} \right)} + \frac{\left. {{s_{I}\Phi_{I}} - {\Phi_{\{{I,M}\}}\left( A_{I} \right)}} \right)^{2}}{2{\sigma_{\{{I,M}\}}^{2}\left( A_{I} \right)}}} \right)}}} \right\rbrack} + \underset{\underset{{epicardium}\mspace{14mu}{shape}\mspace{14mu}{prior}}{︸}}{\int{\int_{\Omega}{{H_{\alpha}\left( \Phi_{O} \right)}\left\lbrack {{\log\left( {\sigma_{\{{O,M}\}}\left( A_{O} \right)} \right)} + \frac{\left( {{s_{O}\Phi_{O}} - {\Phi_{\{{O,M}\}}\left( A_{O} \right)}} \right)^{2}}{2{\sigma_{\{{O,M}\}}^{2}\left( A_{O} \right)}}} \right\rbrack}}}}$where A_(I) is the registration, that is, rigid transformation betweenthe endocardium and its prior shape model [Φ_({I,M}),σ_({I,M})]. Thesame notation can be used for the epicardium.

Integration of visual information modules with anatomical constraintsand prior shape knowledge can be performed using optimization criterioncombining five different components using some weight factors:

${E\left( {\Phi_{I},\Phi_{O},A_{I},A_{O}} \right)} = {\underset{\underset{regularity}{︸}}{w_{1}{E_{G}\left( {\Phi_{I},\Phi_{O}} \right)}} + \underset{\underset{boundary}{︸}}{w_{2}{E_{B}\left( {\Phi_{I},\Phi_{O}} \right)}} + \underset{\underset{region}{︸}}{w_{3}{E_{R}\left( {\Phi_{I},\Phi_{O}} \right)}} + \underset{\underset{anatomy}{︸}}{w_{4}{E_{A}\left( {\Phi_{I},\Phi_{O}} \right)}} + \underset{\underset{{Shape}\mspace{14mu}{Prior}}{︸}}{w_{5}{E_{s}\left( {\Phi_{I},\Phi_{O},A_{I},A_{O}} \right)}}}$The defined objective function refers to two sets of unknown variables.The two sets of unknown variables are two level set functions for theepicardium and endocardium, and two rigid transformations between theevolving representations and the shape prior models. The minimization ofthe objective function can be done using a gradient descend method andthe calculus of variations.

With regard to energy parameters, the boundary term has an insignificantcontribution. Also, the anatomical constraint holds a minor role sinceit varies significantly over space and time. Finally, the regularityconstraint is overwritten by the use of the shape prior term thatprivileges a certain topology. Therefore, two components are used, theregion module and the shape prior term.

The numerical implementation of an embodiment according to the presentinvention is now described. The level set implementation is performedusing a Narrow Band Method. The essence of the Narrow Band Method is toperform level set propagation within a limited zone, that is, aparameter of the DIRAC and HEAVISIDE distributions, located around thelatest position of the propagating contours, in the inward and outwarddirection. Thus, working area is significantly reduced, resulting in asignificant decrease of the computational complexity per iteration.However, the Narrow Band Method requires a frequent re-initialization oflevel set functions that can be performed using a Fast Marchingalgorithm.

An embodiment of the present invention comprises segmentation of theleft ventricle in Magnetic Resonance images by using level setrepresentations.

Another embodiment according to the present invention comprises aformulation that solves segmentation and registration problemssimultaneously. This advantageously allows simultaneously recovering anoptimal segmentation map and a transformation that optimally is a finalresult to a predefined shape model.

Another embodiment according to the present invention, comprises levelset representations, knowledge from medicine and physiology, andavailable visual information for improving the segmentation performanceto physically corrupted data. Furthermore, the mathematical frameworkthat is described, although based on two dimensional visual informationcan be extended to three dimensional space due to the nature of theselected representation.

The teachings of the present disclosure are preferably implemented as acombination of hardware and software. Moreover, the software ispreferably implemented as an application program tangibly embodied on aprogram storage unit. The application program may be uploaded to, andexecuted by, a machine comprising any suitable architecture. Preferably,the machine is implemented on a computer platform having hardware suchas one or more Central Processing Units (“CPUs”), a Random Access Memory(“RAM”), and Input/Output (“I/O”) interfaces. The computer platform mayalso include an operating system and micro instruction code. The variousprocesses and functions described herein may be either part of the microinstruction code or part of the application program, or any combinationthereof, which may be executed by a CPU. In addition, various otherperipheral units may be connected to the computer platform such as anadditional data storage unit and an output unit.

It is to be further understood that, because some of the constituentsystem components and steps depicted in the accompanying drawings may beimplemented in software, the actual connections between the systemcomponents or the process function blocks may differ depending upon themanner in which the present disclosure is programmed. Given theteachings herein, one of ordinary skill in the pertinent art will beable to contemplate these and similar implementations or configurationsof the present disclosure.

Although illustrative embodiments have been described herein withreference to the accompanying drawings, it is to be understood that thepresent disclosure is not limited to those precise embodiments, and thatvarious changes and modifications may be effected therein by one ofordinary skill in the pertinent art without departing from the scope orspirit of the present disclosure. All such changes and modifications areintended to be included within the scope of the present disclosure asset forth in the appended claims.

1. A method for segmentation of medical images comprising the steps of:inputting a mathematical functional into a system for tracking movinginterfaces, wherein said mathematical functional accounts forglobal/local shape properties of a physiological object of interest thatis to be recovered; visualizing image information wherein saidinformation includes boundary and regional data; and recovering saidphysiological object of interest by optimizing said mathematicalfunctional integrated with said image information.
 2. The method ofclaim 1, further comprising the step of combining said mathematicalfunctional with a level set objective function having a medical shapemodel.
 3. The method of claim 2, wherein said level set objectivefunction includes a bi-directional boundary flow.
 4. The method of claim2, wherein said medical shape model includes a medical shape modelhaving a degree of variability.
 5. The method of claim 1, furthercomprising the step of combining said mathematical functional with aphysiology based functional for constraining the solution space and aterm that accounts for shape driven consistency.
 6. The method of claim1, further comprising the step of combining said mathematical functionalwith an intensity regional function that maximizes a posteriorisegmentation probability.
 7. The method of claim 1, further comprisingthe step of optimizing a resulting objective function using a gradientdescent method.
 8. A method for segmentation of medical imagescomprising the steps of: developing a medical shape model on level setspace for tracking moving interfaces; using said medical shape model forintroducing a physiological shape prior in an energetic form; andrecovering a physiological object of interest by minimizingnon-stationary distance between an evolving interface and said medicalshape model.
 9. The method of claim 8, wherein said medical shape modelis developed using a variational framework to create a non-stationarypixel-wise model that accounts for shape variabilities.
 10. The methodof claim 8, further comprising the step of minimizing non-stationarydistance between an evolving interface and said medical shape model. 11.The method of claim 8, further comprising the step of integrating saidmedical shape model and said physiological shape prior in energetic forminto a data-driven variational method that performs image segmentation.12. The method of claim 8, further comprising the step of integrating anintensity-based regional module that maximizes a posteriori segmentationprobability.
 13. The method of claim 8, further comprising the step ofintegrating a physiology-based module that constrains solution space anda term that accounts for shape driven consistency.
 14. The method ofclaim 8, further comprising the step of optimizing a resulting objectivefunction using a gradient descent method.
 15. A program storage devicereadable by machine, tangibly embodying a program of instructionsexecutable by the machine to perform method steps for segmentation ofmedical images, the method steps comprising: developing a medical shapemodel on level set space for tracking moving interfaces; using saidmedical shape model for introducing a physiological shape prior in anenergetic form; and recovering a physiological object of interest byminimizing non-stationary distance between an evolving interface andsaid medical shape model.
 16. The program storage device of claim 15,wherein said medical shape model is developed using a variationalframework to create a non-stationary pixel-wise model that accounts forshape variabilities.
 17. The program storage device of claim 15, whereinthe method steps further comprise the step of minimizing non-stationarydistance between an evolving interface and said medical shape model. 18.The program storage device of claim 15, wherein the method steps furthercomprise the step of integrating said medical shape model and saidphysiological shape prior in energetic form into a data-drivenvariational method that performs image segmentation.
 19. The programstorage device of claim 15, wherein the method steps further comprisethe step of integrating an intensity-based regional module thatmaximizes a posteriori segmentation probability.
 20. The program storagedevice of claim 15, wherein the method steps further comprise the stepof integrating a physiology-based module that constrains solution spaceand a term that accounts for shape driven consistency.
 21. The programstorage device of claim 15, wherein the method steps further comprisethe step of optimizing a resulting objective function using a gradientdescent method.